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We investigate numerically the failure process when two 
elastic media, one hard and one soft that have been glued 
together thus forming a common interface, are pulled apart. 
We present three main results: (1) The area distribution of 
simultaneously failing glue (bursts) follows a power law con- 
sistent with the theoretically expected exponent 2.5, (2) the 
maximum load and displacement before catastrophic failure 
scale as and respectively, where L is the linear size of 
the system, and (3) the area distribution of failed glue regions 
(clusters) is a power law with exponent —1.6 when the system 
fails catstrophically. 

PACS number(s): 83.80.Ab, 62.20.Mk, 81.40. Np 



I. INTRODUCTION 

The failure of interfaces under stress is a problem that 
has obvious important technological relevance. In addi- 
tion, from a more fundamental point of view, this prob- 
lem exhibits very interesting features. It is the aim of this 
paper to bring out some of these features by means of a 
numerical model based on a discretization of the original 
problem. 

During more than a decade, failure processes in dif- 
ferent contexts have caught the attention of the physics 
community. For a considerable longer period, the me- 
chanics community has been involved in the study of 
such phenomena. In order to place the present study 
in its proper context, we need to go back to 1926 with 
the study of Peirce on what today is known as the global 
load sharing fiber bundle |]l| . This consists of N parallel 
fibers, each with its own breaking threshold and con- 
nected in such a way that when a fiber fails, the load it 
was carrying would be dsitributed equally among all the 
surviving fibers. In 1945, Daniels published a very thor- 
ough study of this model, which today forms the starting 
point of any excursion into this field |^ . The model has 
since these early days been generalized in many direc- 
tions, one of which consists in replacing the "democratic" 
load-sharing rule by different local ones. One much stud- 
ied variant is the local load sharing model, where the load 



on the failing fiber is distributed equally among the near- 
est surviving fibers The average load-deformation 
characteristics of the global load sharing model was cal- 
culated by Daniels 1^. Corresponding work on the local 
load sharing model may be found in Refs. 1^-^. There 
has also been several studies of time-dependent phenom- 
ena in connection with the two variants of the fiber bun- 
dle model, see Ref. §. (This paper in addition contains 
a very thorough review of the literature in this field.) 
There are also a number of studies "on the market" that 
may be placed between the two extremes that the global 
and local load sharing models consitute. Among them, 
we find the early study by Newman and Gabrielov [Q, 
who constructed a hierarchically connected fiber bundle. 
Other work on hierarchical fiber bundle models may be 
found in Refs. | -|l^. 

Much work by the physics community has gone into 
studying network models, of which the fuse model is 
the most well known ||ll|,|l^. This model consists of a 
network of electrical fuses where their burn-out thresh- 
olds have been drawn from some probability distribu- 
tion. This model may be regarded as yet another gen- 
eralization of the fiber bundle model, however, this time 
along the axis on which we find chains of fiber bundles 
|l3| . Among the several interesting questions that have 
been studied in connection with the fuse model, we men- 
tion the question of whether the breakdown process has 
the character of a second or first order phase transition 
|l^-^. Central to this question is the question of the 
distribution of fuses that burn out simultaneously or — 
equivalently in the fiber bundle model — the number 
of fibers that fail simultaneously. This question was first 
raised and solved analytically in the context of the global 
load sharing fiber bundle ||l^ and then for the local load 
sharing model . The same question was first studied 
in connection with the fuse model in Ref. . 

The particular problem we study here, elastic interfa- 
cial failure, has been addressed in the literature earlier 
by Delaplace et al. p0| , pl| , p2[ . The system consists of two 
elastic media that have been welded together, thus shar- 
ing a common interface. In general, the media can have 
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different elastic constants. However, for the sake of sim- 
plicity and without loss of generality, we assume one of 
the media to be infinitely stiff while the other is elastic. 
We can view this simplification as an effective represen- 
tation of the original system since it does not change the 
physics. Furthermore, we assume that the "soft" medium 
is uniform with respect to its elastic properties; it has the 
same elastic properties everywhere. However, the local 
strength of the glue — defined as the maximum local 
load it may sustain without failing — varies from point 
to point along the interface. This is the source of disorder 
in the system. In real systems this disorder would typi- 
cally be correlated. In this first attack on the problem, 
we assume the disorder to be uncorrelated. Our main 
interest here is to understand how correlations develop 
due to the failure process. In Section || we describe the 
numerical model in detail. 

The two joined media are subjected to a progressive 
uniform load perpendicular to the glue interface. Lo- 
cal failures will develop in the interface which changes 
the stress field on the remaining intact interface. These 
changes in the stress field will compete with the local 
strength of the glue to determine where the next fail- 
ure happens. Sometimes, a local failure will occur due 
to the glue being particularly weak at that point on the 
interface, other times failure will occur due to enhance- 
ments in the local stress field. This competition leads to 
the development of spatial correlations both in the stress 
field and in the failure patterns, and in the failure process 
itself. 

The two media can be pulled apart by controlling (fix- 
ing) either the applied force or the displacement. The 
displacement is defined as the change in the distance be- 
tween two points, one in each medium, positioned far 
from the glued interface. Clearly, the line connecting 
these points is perpendicular to the average position of 
the interface. In our case, the pulling is accomplished 
by controlling the displacement. As the displacement is 
increased very slowly, glued points will fail. Sometimes 
the failed regions are very small, other times, the failed 
region is larger. Such events, when a large area fails "in- 
tantaneously" compared to the time scale at which the 
displacement is changed, are called bursts. One of the 
quantities of interest to us is the burst distributio n fl^] 
as the failure process evolves. We show in Section [II A 
that this distribution follows a power law. 



In Section III E , we investigate the scaling properties 
of the load and displacement of the system at the point 
when the failure process becomes unstable. This is the 
point at which any further increase of either load or dis- 
placement will lead to a catastrophic burst where all re- 
maining glue fails. This point defines the strength of the 
interface, and the question we pose is how this scales with 
the system size. 

We then investigate the geometrical properties of the 
failed regions at the point when catastrophic failure oc- 



curs m Section |IIIC| . We find that the area distribution 
of the failed regions follow a power law. 

We present our conclusions and outlook for further 
work in Section |v| . 



II. MODEL 

The system described in the Introduction is continu- 
ous: Two media, one elastic, the other infinitely stiff, 
are glued together thus forming a common interface. In 
order to treat this problem numerically, the continuum 
problem is replaced by a discrete one. We use a discrete 
model for this. In Section II A we describe the discrete 



model, while the numerical algorithms are discussed in 
Section flB. 



A. Description of Model 

We discretize the glued interface by replacing it with 
two two-dimensional square L x L lattices with periodic 
boundary conditions. The lower one represents the hard, 
stiff surface and the upper one the elastic surface. The 
nodes of the two lattices are matched (i.e. there is no 
relative lateral displacement). The glue is modelled by 
springs connecting opposing nodes in the two lattices. 
These harmonic springs all have the same spring constant 
(set to unity) but breaking thresholds randomly drawn 
from a uniform distribution between zero and one. As 
the two glued media are separated, the forces carried by 
the springs increases. When the force carried by a spring 
reaches its breaking threshold, it breaks irreversibly and 
the forces redistribute. The springs are thus broken one 
by one until the two media are no longer in mechanical 
contact. As this process is proceeding, the elastic body 
is of course deforming to accomodate the changes in the 
forces acting on it. 

The equations governing the system are as follows. The 
force, fi, carried by the zth spring is given by Hooke's law, 



f^ 



(1) 



where k is the spring constant, D is the displacement of 
the hard medium, and Ui is the deformation of the elastic 
medium at site i. All unbroken springs have k — 1 while 
a broken spring has fc = 0. The quantity (ui — D) is, 
therefore, the length spring i was stretched. In addition, 
a force applied at a point on an elastic surface will de- 
form this surface over a region whose extent depends on 
its elastic properties. This is described by the coupled 
system of equations. 



(2) 



where the elastic Green function, is given by p3 24 
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_ 1^ dx' dy' 

nEa^ y„,/2 \{x-x',y-y')\ ' 

In this equation, s is the Poisson ratio, e the elastic con- 
stant, and |j — j| the distance between sites i and j. The 
indices i and j run over all sites. Strictly speaking, 
this Green function applies for a medium occupying the 
infinite half space. However, with a judicious choice of 
elastic constants, we may use it for a finite medium if its 
range is small compared to L, the size of the system. 
By combining equations (]l|) and (|^), we obtain 

(I + KG)/ = Ki3, (4) 

where we are using matrix-vector notation. I is the 
X identity matrix, and G is the Green function 
represented as an x dense matrix. The constant 
vector D is dimensional. The diagonal matrix K is 
also X L^. Its matrix elements are either 1, for unbro- 
ken springs, or for broken ones. Of course K and G do 
not commute (except initially when there are no broken 
springs). 

Once equation (^) is solved for the force, /, equation 
easily yields the deformations of the elastic surface. 

B. Numerical Method: Fourier Acceleration 

Equation is of the familiar form Ax = b. Since the 
Green function connects all nodes to all other nodes, the 

X matrix A is dense which puts severe limits on the 
size of the system that may be studied. There are direct, 
time consuming methods to deal with such matrices, see 
Ref. However, as we shall see, this problem may 

be circumvented and much more efficient methods may 
be employed such as the Conjugate Gradient algorithm 

(CG) 

The simulation proceeds as follows: We start with all 
springs present, each with its randomly drawn breakdown 
threshold. The two media are then pulled apart, the 
forces calculated using CG, and the spring which is the 
nearest to its threshold is broken, i.e. the matrix element 
corresponding to it in the matrix K is zeroed. Then the 
new forces are calculated, a new spring broken and so on 
till all springs have been broken and the media separated. 

However, there are two problems that render the sim- 
ulation of large systems extremely difficult. The first is 
that since G is x dense matrix, the number of oper- 
ations per CG iteration scales like L'^. Even more serious 
is the fact that as the system evolves and springs are bro- 
ken, the matrix (I -I- fcG) becomes very ill-conditioned. 

To overcome the problematic L"* scaling of the algo- 
rithm we note that the Green function is diagonal in 
Fourier space. Consequently, doing matrix-vector multi- 
plications using FFTs the scaling is much improved and 
goes like ln(L). Symbolically, this can be expressed as 
follow: 




color scheme indicates when in the failure process a given 
bond failed, the lighter, the earlier. The lattice was 128 x 128 
with an elastic constant e = 10. 

(I + KF"^FG)F"^F/ = KZ? , (5) 

where F is the EFT operator and F^-*- its inverse 
(F^-'-F = 1). Since I and K are diagonal, operations 
involving them are performed in real space. With this 
formulation, the number of operations/iteration in the 
CG algorithm now scales like ln(L). 

To overcome the runaway behavior due to the ill- 
conditioning we need to precondition the matrix [ p6| . 
This means that instead of solving equation (j^), we solve 
the equivalent problem 

Q(I + KF-^FG)F-^F/ = QKl3 , (6) 

where we simply multiplied both sides by the arbitrary, 
positive definite preconditioning matrix Q. Clearly, the 
ideal choice is Qo = (I+KG)^^ which would always solve 
the problem in one iteration. Since this is not possible 
in general, we look for a form for Q which satisfies the 
following two conditions: (1) As close as possible to Qo, 
and (2) fast to calculate. The choice of a good Q is 
further complicated by the fact that as the system evolves 
and springs are broken, corresponding matrix elements of 
K are set to zero. So, the matrix (I -I- KG) evolves from 
the initial form (I -t- G) to the final one I. We were 
not able to find a fixed Q that worked throughout the 
breaking process. 

We therefore chose the form 

Q = I - (KG) + (KG)(KG) - (KG)(KG)(KG) + ... 

(7) 

which is nothing but the Taylor series expansion of 
Qo = (I + KG)^^. For best performance, the number of 
terms kept in the expansion is left as a parameter since 
it depends on the physical parameters of the system. It 
is important to emphasize the following points, (a) As 
springs are broken, the preconditioning matrix evolves 
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2. Distance between succesively broken bonds. The 
was 128 X 128 with an elastic constant e = 10. 



with the ill-conditioned matrix and, therefore, remains a 
good approximation of its inverse throughout the break- 
ing process, (b) All matrix multiplications involving G 
are done using FFTs. (c) The calculation of Q can be 
easily organized so that it scales like nL^ ln(i) where n is 
the number of term kept in the Taylor expansion, equa- 
tion (|). 

We therefore have a stable accelerated algorithm which 
scales essentially as the volume of the system. For exam- 
ple, for a 128 x 128 system, and taking n — 5, the CG 
algorithm always converged in four or five iterations with 
the prescribed precision of 10^^-^. 



III. RESULTS 
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FIG. 3. Distance between succesively broken bonds. The 
lattice was 128 x 128 with an elastic constant e = 100. 
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FIG. 4. Same as in Fig. nl but with a narrow uniform 
threshold distribution in the interval [9.5, 10.5]. 



We now present the results of our numerical simula- 
tions. We show in Fig. |^ a representation of the failure 
process. Each elementary square represents a spring (a 
bond), and the gray scale indicates when a particular 
spring failed: The darker the color, the earlier the fail- 
ure. In this particular case, the elastic constant e was set 
to 10. There are no apparent spatial correlations between 
the failing bonds in this figure. However, we show in Fig. 
^ the distance between successively failing bonds for the 
same disorder realization of Fig. ^ We see clearly in this 
figure that early in the process there is no localization 
effect: Bonds tend to break far apart, the location being 
determined by the strength of bonds, i.e. early failure 
is disorder driven. However, halfway into the breakdown 
process, localization clearly develops. In Fig. ||, we show 
the corresponding plot for e = 100. In this case localiza- 
tion never develops for this size system and destribution 
of thresholds. 



On the other hand, if the threshold distribution is 
much narrower than [0, 1] used above, localization can de- 
velop early. For example, we show in Fig. ^ the fracture 
graymap (like Fig. |^) for a uniform threshold distribu- 
tion in the interval [9.5, 10.5]. We clearly see the fracture 
starting towards the center of the figure and spreading 
out in a spiral till finally the symmetry is broken and the 
system ruptures along one of the lattice directions. 

Fig. H shows the force-displacement curve for a system 
with elastic constant e = 10. Whether we control the 
applied force, or the displacement, -D, the system will 
eventually suffer catastrophic collapse. However, this is 
not so when e = 100 as shown in Fig. ^. In this case, only 
controlling force will lead to catastrophic failure. In the 
limit when e ^ cx), the model becomes the global load- 
sharing fiber bundle model ||l|,||, where F = {1 — D)D. In 
this limit there are no spatial correlations and the force 
instability is due to the the decreasing total elastic con- 
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FIG. 5. Force-displacement curve, 128 x 128 systems with 
e = 10. 
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FIG. 6. Force-displacement curve, 128 x 128 systems with 
e = 100. 



stant of the system making the force on each surviving 
bond increase faster than the typical spread of threshold 
values. No such effect exists when controlling displace- 
ment D. However, when the elastic constant, e, is small, 
spatial correlations in the form of localization do develop, 
and these are responsible for the diplacement instability 
which is seen in Fig. H. In other words, the localization 
clearly visible in Fig. g starts to develop when the system 
is near the peak of its force-displacement curve, and dom- 
inates when the system is on the negative slope branch 
of that curve. 



A. Burst Distribution 

We now turn to the study of the burst distribution. We 
define the size of a burst. A, in our model as the number 
of bonds that fail simultaneously while the force F is held 
constant. In the global load-sharing fiber bundle model 
it has been shown that the burst distribution is given by 
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7V(A,i5) = — n(A"(x-Xe)) 



(8) 



where Xc is the damage, i.e. the density of broken bonds, 
at which the model fails catastrophically, n is a crossover 
function which approaches a constant when the argument 
approaches zero, and which falls off as exp(— y^) as the 
argument y is large. Furthermore, 



and 



1 

(T = — 

2 



(9) 



independent of the threshold distribution. 

We show in Figs. |^ and ^ the burst distribution for 
e = 10 and 100. In both cases we find that the burst 



distribution follows a power law with an exponent t — 
—2.6 ± 0.1. We may argue that the exponent is the same 
as the one found in the global loading fiber bundle, Eq. 
®; T = 5/2 in the following way. The characteristics, 
F = F{x) must have a quadratic maximum somewhere. 
For e — 100, such a maximum exists in the middle of 
failure process as seen in Fig. ||, whereas for e — 10, 
the system only approaches such a maximum near to- 
tal failure, see Fig. |^. Assuming that the fluctuations 
about the average characteristics are brownian — which 
can be shown analytically in the limit e ^ oo p7|-p9|] 
— near the maximum the probability to find a burst of 
size A is proportional to A~^/^ exp[— A(a: — Xc)^]. This 
result comes from a mapping onto the Gambler's ruin 
problem [^ . Furthermore, in order to guarantee that 
the burst is not a burst within an even larger burst, the 
starting point of the burst must be the highest point on 
the characteristics that has occured so far in the fail- 
ure process. This condition may also be mapped onto 
the Gambler's ruin problem, and leads to an extra fac- 
tor {x — Xc) in the probability for a burst to occur. The 
probability to find a burst of size A throughout the fail- 
ure process is then the integral over x as a; approaches 
Xc, Z^'' dx{x — a;c)A~'^/^ exp[— A(x — Xc)"^] which is pro- 
portional to A~^/^. 

As can be seen in Figs. ^ and ||, the numerical data are 
consistent with the expected value r = 5/2. 



B. Strength Scaling 

The load-displacement curves for different system sizes 
L coincide when we use the reduced variables F/ L° and 
D/V' where a — 2.0 and h = 0.0, as seen in Fig. ||. 
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FIG. 7. Burst distribution for 128 x 128, e 
of the straight Une is —2.5. 
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FIG. 9. Scahng of the failure characteristics for systems 
with L = 128(e = 10),64(e = 5) and 32(e = 2.5) using the 
reduced variables F/L^ ° and D/L° °. 

We expect the exponent a = 2 since F/L^ is the normal 
stress on the surface. In the case of an infinitely stiff 
system, we expect 5 = 0. The elastic system studied here 
behaves in the same way as long as the elastic constant, 
e, is also scaled with L. For example for L = 128 we took 
e — 10, for L = 64 e should take half that value in order 
to reproduce the physics. This is easy to understand 
considering the dependence of the Green fucntion, Eq. |^, 
on the elastic constant. 




Burst size 

FIG. 8. Burst distribution for 128 x 128, e = 100. The 
slope of the straight line is —2.5. 



C. Spatial Damage Distribution at Failure 

As the failure process proceeds, there is an increasing 
competition between local failure due to stress enhance- 
ment and local failure due to local weakness of material. 
As we saw above, when we control the displacement, D, 
and e is sufficiently small (for example e = 10), catas- 
trophic failure eventually occurs due to localization. The 
onset of this localization, i.e. the catastrophic regime, 
occurs when the two mechanisms are equally important. 
One may suspect that criticality due to self organiza- 
tion 1^ occurs at this point. In order to test whether 
this is the case, we have measured the size distribution 
of broken bond clusters at the point when D reaches its 
maximum point on the F — D characteristics, i.e. the 
onset of localization and catastrophic failure. The anal- 
ysis was performed using a Hoshen-Kopelman algorithm 
We show the result in Fig. |o|, for 56 disorder re- 
alizations, L = 128 and e = 10. The result is consistent 
with a power law distribution with exponent —1.6, and 
consequently with self organization. 



6 



10^ 



10' 



0) 

o 

0) 



f 10" 



10-^ 





- 


- 

r 128X128 


: 


e=10 


"••vv: 


slope = -1.6 




56 realizations 










10"' 

10" 10' 10' 10" 

Cluster Size 

FIG. 10. Area distribution of zones where glue has failed 
for systems of size 128 x 128 and elastic constant e = 10. The 
straight line is a least square fit and indicates a power law 
with exponent —1.6. e = 10 



IV. CONCLUSION 

We have studied numerically the failure of the glued 
interface between an elastic and an infinitely stiff blocks 
of material. To this end we have developed a new, stable 
and accelerated algorithm which scales L^ln(I/) which 
enabled us to study much bigger systems than previously 
possible. 

Our main physical results arc: (1) the distribution of 
simultaneously failing glue (bursts) is a power law with 
exponent —2.6 ± 0.1, which is consistent with the burst 
distribution found in the global load sharing fiber bun- 
dle problem. (2) The point of catastrophic failure scales 
as in force and in displacement. (3) The area 
distribution of failed regions (clusters) at the onset of 
catastrophic failure when displacement is the control pa- 
rameter is consistent with a power law with an exponent 
equal to —1.6. This hints at self organization. 

In addition, we saw that for large e, e.g. e = 100, the 
system does not suffer catastrophic failure, and there is 
no localization. On the other haand, smaller values of 
e, e.g. e = 10, resulted in catastrophic failure due to 
localization. By doing the simulations for various values 
of e we estimate that failure due to localization starts 
to occur for e ~ 35 — 40. As we will see below, these 
values of e obtained for 128 x 128 systems should be scaled 
appropriately when the size of the system is changed. 

The disorder in our system was uncorrelated. As men- 
tioned above, it is realistic to introduce correlations as 
exist for example in fracture surfaces. This can be done 
by generating spring breaking thresholds that have the 
desired correlations. Furthermore, we have used a flat 



distribution for the disorder. One can use other distri- 
butions, e.g. r", where r is a uniformly distributed ran- 
dom number and a an exponent that can be negative. It 
is known from random fuse models of fracture that the 
breakdown process depends on the value of a. It is not 
clear how these issues will modify our current results. 
This work is in progress. 

Another work in progress is to study the propagating 
fracture front when the above glued media are ripped 
apart by pulling only on one side of the L x L square sys- 
tem. As before the breaking thresholds can be correlated 
or uncorrelated. These results will also be compared with 
the results of experiment currently underway. 

Finally, we have chosen to introduce disorder into the 
breaking thresholds of the springs. However, we can just 
as easily introduce it into the spring constants them- 
selves, again with or without correlations. 
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